//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry                         : env klmperry

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"

// base
cd  $perrydatas
use perry_base.dta, clear
gen     treat = 1 if C2 == 2
replace treat = 0 if C2 == 1
gen perry = 1
keep    id treat sb_iq_5 perry
rename sb_iq_5 earlyiq
tempfile perry
save "`perry'", replace

// abc
cd $abcjpeanalysis
use append-abccare_iv.dta, clear
keep if program == "abc"
egen sb45y = rowmean(sb4y sb5y)
rename sb45y earlyiq
keep id treat earlyiq
gen abc = 1
tempfile abc
save "`abc'", replace

// ihdp
cd $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs
rename ihdp id
rename tg treat
rename iqcage earlyiq
gen    ihdp = 1 
keep   id treat earlyiq ihdp twin bwg site 
append using "`perry'"
append using "`abc'"
gen control = 1 - treat

// specs for run
gen ihdps = 1 if ihdp == 1 & twin == 0
gen ihdpt = 1 if ihdp == 1 & twin == 1

// strata
global perry_strata
global abc_strata
global ihdps_strata strata(bwg site treat)
global ihdpt_strata strata(bwg site treat)

// treatment effects
// programs
foreach prog in abc perry ihdps ihdpt {
	bootstrap, ${`prog'_strata} reps($bootstraps): reg earlyiq treat if `prog' == 1
	matrix b = e(b)
	matrix V = e(V)
	matrix se`prog' = sqrt(V[1,1])
	matrix md`prog' = b[1,1]
	matrix  t`prog' = abs(md`prog'[1,1]/se`prog'[1,1])
	matrix df`prog' = e(N) - e(rank)
	matrix  p`prog' = 2*(1 - normal(t`prog'[1,1]))
	matrix  n`prog' = e(N)
	
	reg earlyiq control treat if `prog' == 1, nocons
	matrix b = e(b)
	matrix bc = b[1,1]
	matrix bt = b[1,2]
	
	matrix `prog' = [bc,bt,md`prog',p`prog',n`prog']
}
matrix all_0 = [perry \ abc \ ihdps]
matrix all_1 = [perry \ abc \ ihdpt]

// ihdp by site
keep if ihdp == 1
// sites 

gen      states = 1 if site == "ARK"
replace  states = 5 if site == "EIN"
replace  states = 4 if site == "HAR"
replace  states = 3 if site == "MIA"
replace  states = 6 if site == "PEN"
replace  states = 7 if site == "TEX"
replace  states = 8 if site == "WAS"
replace  states = 2 if site == "YAL"

levelsof states, local(sites)
foreach sib of numlist 0 1 {
	foreach state in `sites' {
		bootstrap, strata(bwg site treat) reps($bootstraps): reg earlyiq treat if states == `state' & twin == `sib'
		matrix b = e(b)
		matrix V = e(V)
		matrix se = sqrt(V[1,1])
		matrix md = b[1,1]
		matrix  t = abs(md[1,1]/se[1,1])
		matrix df = e(N) - e(rank)
		matrix  p = 2*(1 - normal(t[1,1]))
		matrix  n = e(N)
		
		reg earlyiq control treat if states == `state' & twin == `sib', nocons
		matrix b = e(b)
		matrix bc = b[1,1]
		matrix bt = b[1,2]
		
		matrix all_`sib' = [all_`sib' \ [bc,bt,md,p,n]]
	}
}

matrix colnames all_0 = bc_0 bt_0 md_0 p_0 n_0
matrix colnames all_1 = bc_1 bt_1 md_1 p_1 n_1
matrix all = [all_0,all_1]

clear
svmat all, names(col)
// labels

foreach sib of numlist 0 1 {
	tostring n_`sib', replace force format(%15.0fc)
	replace  n_`sib' = "N = " + n_`sib'
	tostring md_`sib', gen(md_`sib'str) force format(%15.1fc)
	replace  md_`sib'str = "{&Delta} = " + md_`sib'str
}
	

gen n = _n*2 - 2
gen ppt1 = 116
gen ppt2 = 114

// plots
global la_0 Singletons
global la_1 Twins

foreach num of numlist 0 1 {
	#delimit
	twoway     (bar     bc_`num'  n, color(gs10) barwidth(.8)              )
		   (bar     bt_`num'  n, fcolor(none) lcolor(gs0) barwidth(.8) xline(5, lwidth(vthick) lcolor(gs10)))

		   (scatter ppt1  n, msymbol(none) mlabel(md_`num'str)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
		   (scatter ppt2  n, msymbol(none) mlabel(n_`num')        mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
		   
		   (scatter bt_`num'    n    if  p_`num' < 0.05                  , msize(small)  mcolor(gs0) msymbol(circle))
		   (scatter bt_`num'    n    if  p_`num' >= .05 & p_`num'    <= 0.10 , msize(small)  mcolor(gs0) msymbol(square))
		  

		   
		   ,
			  legend(label(1 "Control") label(2 "Treatment")
				 label(5 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
				 label(6 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
				cols(2) rows(2) order(1 2 5 6) size(medsmall) position(6) region(lcolor(gs0)))
					 
			  ylabel(70[12]118, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
			  xlabel(0 "{bf:Perry}" 2 "{bf:ABC}" 4 `" "${la_`num'}" "{bf:IHDP}" "'
				 6  "AR" 8  "CT" 10 "FL" 12 "MA" 13 `" " " "{bf:IHDP}" "'
				 14 "NY" 16 "PA" 18 "TX" 20 "WA"
			  ,  labsize(small) grid glcolor(gs14) glpattern(solid))
			  ytitle("Unadjusted Mean", size(small))
			  xtitle("")
			  graphregion(color(white)) plotregion(fcolor(white)); 
	#delimit cr
	cd $output
	graph export TE_sbiq_by_site_s`num'.eps, replace
}
